R Notebook
Machine learning project
Authors: Jose Pérez Cano & Álvaro Ribot Barrado
0. Libraries
install.packages("klaR")
install.packages("TunePareto")
install.packages("rgl")
install.packages("glmnet")
install.packages("ca")# LDA/ QDA
library(MASS)
# RDA
library(klaR)
# Multinomial
library(nnet)
# Cross-Validation
library(TunePareto)
# Naive Bayes
library(e1071)
# k-NN
library(class)
# 3d
library(rgl)
# LASSO
library(Matrix)
library(glmnet)Loading required package: foreach
Loaded glmnet 2.0-16# Correspondence analysis
library(ca)1. Read data
set.seed(2105)
setwd("../data")The working directory was changed to /Users/joseperezcano/Desktop/CFIS segundo curso/AA1/Project/data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.clev <- read.csv("cleveland.csv", header=F)
hung <- read.csv("hungarian.csv", header=F)
va <- read.csv("long-beach-va.csv", header=F)
switz <- read.csv("switzerland.csv", header=F)
clev$location <- "cleveland"
hung$location <- "hungarian"
va$location <- "long-beach-va"
switz$location <- "switzerland"
heart1 <- rbind(clev, hung)
heart2 <- rbind(va, switz)
heart <- rbind(heart1, heart2)
head(heart)2. Preprocess data
We apply clustering and several plotting techniques to have an idea of the dataset. In case there are NAs we will use k-NN for imputation.
To extract new features we will apply PCA and FDA and keep this new features and components apart.
# It says which columns are all missings
# The index are returned in negative to eliminate them
na.columns <- function(dd){
rmlist <- c()
for (i in 1:ncol(dd)){
if (min(dd[,i]) == max(dd[,i]) & min(dd[,i])==-9){
rmlist <- c(rmlist, i)
}
}
-rmlist
}
clev <- clev[,na.columns(clev)]
# Returns columns with more NA than a given threshold, also in negative
much.na.cols <- function(dd, threshold){
rmlist <- c()
for (i in 1:ncol(dd)){
if (sum(dd[,i]==-9) > threshold){
rmlist <- c(rmlist, i)
}
}
-rmlist
}
clev <- clev[, much.na.cols(clev, 60)]
# Applies k-nearest neighbour imputation for a given variable
knn.imputation = function (dd, variable, varname, k)
{
aux = subset (dd, select = names(dd)[names(dd) != varname])
aux1 = aux[!is.na(variable),]
aux2 = aux[is.na(variable),]
# Neither of aux1, aux2 can contain NAs
knn.inc = knn (aux1,aux2, variable[!is.na(variable)], k)
variable[is.na(variable)] = knn.inc
variable
}
# This are the variables which values where substituted by dummy values.
dummy <- c("V1", "V2", "V36", "V69", "V70", "V71", "V72", "V73", "V28", "location")
clev <- clev[,!(names(clev) %in% dummy)]
# knn imputation for clev
na.names <- names(clev)[-much.na.cols(clev, 0)]
for (name in na.names){
clev[, name][clev[, name] == -9] <- NA
clev[, name] <- knn.imputation(clev, clev[,name], name, 7)
}Now, we analyse the correlations among variables since it will have impact on later models.
corr.factors <- cor(clev)
which(abs(corr.factors)-diag(diag(corr.factors))>0.9, arr.ind=T) row col
V55 34 12
V57 36 14
V31 21 20
V29 20 21
V20 12 34
V22 14 36rm.correlated <- c("V57", "V55")
clev <- clev[,!(names(clev) %in% rm.correlated)]
factores <- c("V58", "V4", "V9", "V16", "V18", "V19", "V20", "V21", "V22", "V23", "V24", "V25", "V26", "V27", "V38", "V39", "V41", "V51", "V56", "V11", "V59", "V60", "V61", "V63", "V65", "V67", "V68")
for (f in factores){
clev[,f] <- as.factor(clev[,f])
}
# Dummy level gets replaced
clev$V25[clev$V25 == 2] <- 1
clev$V25 <- droplevels(clev$V25)
levels(clev$V25)[1] "0" "1"This is the dataset after the treatment for missings.
summary(clev)
V3 V4 V9 V10 V11
Min. :29.00 0: 91 1: 22 Min. : 94.0 0:108
1st Qu.:48.00 1:191 2: 43 1st Qu.:120.0 1:174
Median :55.00 3: 84 Median :130.0
Mean :54.41 4:133 Mean :131.6
3rd Qu.:61.00 3rd Qu.:140.0
Max. :77.00 Max. :200.0
V12 V14 V15 V16 V18
Min. :126.0 Min. : 0.00 Min. : 0.00 0:240 0:107
1st Qu.:213.0 1st Qu.: 0.00 1st Qu.: 0.00 1: 42 1:175
Median :244.0 Median :10.00 Median :15.00
Mean :249.1 Mean :16.64 Mean :15.08
3rd Qu.:277.0 3rd Qu.:30.00 3rd Qu.:30.00
Max. :564.0 Max. :99.00 Max. :54.00
V19 V20 V21 V22 V23 V24
0:138 1 :38 4 : 13 81:67 0:271 0:186
1: 2 12 :33 15 : 13 82:96 1: 11 1: 95
2:142 2 :31 20 : 13 83:86 2: 1
8 :30 13 : 12 84:33
6 :26 16 : 12
11 :25 21 : 12
(Other):99 (Other):207
V25 V26 V27 V29 V31
0:211 0:252 0:248 Min. : 1.800 Min. : 3.000
1: 71 1: 30 1: 34 1st Qu.: 6.500 1st Qu.: 7.000
Median : 8.500 Median : 9.000
Mean : 8.418 Mean : 9.754
3rd Qu.:10.075 3rd Qu.:12.000
Max. :15.000 Max. :18.000
V32 V33 V34
Min. : 71.0 Min. : 40.00 Min. : 84.0
1st Qu.:133.2 1st Qu.: 65.00 1st Qu.:154.0
Median :153.5 Median : 74.00 Median :168.0
Mean :149.8 Mean : 75.12 Mean :168.1
3rd Qu.:165.8 3rd Qu.: 84.00 3rd Qu.:184.0
Max. :202.0 Max. :119.00 Max. :232.0
V35 V37 V38 V39
Min. : 26.00 Min. : 50.00 0:190 0:276
1st Qu.: 70.00 1st Qu.: 80.00 1: 92 1: 6
Median : 80.00 Median : 85.00
Mean : 78.74 Mean : 84.95
3rd Qu.: 85.00 3rd Qu.: 90.00
Max. :120.00 Max. :110.00
V40 V41 V43 V44 V51
Min. :0.000 1:135 Min. : 24.0 Min. :0.0000 1: 2
1st Qu.:0.000 2:129 1st Qu.: 92.0 1st Qu.:0.0000 3:159
Median :0.800 3: 18 Median :118.0 Median :0.0000 6: 14
Mean :1.027 Mean :123.6 Mean :0.6702 7:107
3rd Qu.:1.600 3rd Qu.:152.8 3rd Qu.:1.0000
Max. :6.200 Max. :270.0 Max. :3.0000
V56 V58 V59 V60 V61 V63 V65
5 : 14 0:157 1:270 1:242 1:224 1:238 1:236
21 : 14 1: 50 2: 12 2: 40 2: 58 2: 44 2: 46
14 : 12 2: 31
1 : 11 3: 32
17 : 11 4: 12
30 : 11
(Other):209
V67 V68
1:233 1:246
2: 49 2: 36
2.1 Visualizations
par(mfrow = c(2,3))
for(i in 1:ncol(clev)){
if (!is.factor(clev[,i])) {
hist(clev[,i], main = names(clev)[i], xlab="Values")
}
}
par(mfrow = c(2,3))
for(i in 1:ncol(clev)){
if (!is.factor(clev[,i])) {
boxplot(clev[,i], xlab = names(clev)[i])
}
}
par(mfrow = c(3,3))
for(i in 1:ncol(clev)){
if (is.factor(clev[,i])) {
hist(as.numeric(as.character(clev[,i])), main = names(clev)[i], xlab="Values")
}
}
names_num <- c()
for(i in 1:ncol(clev)){
if (!is.factor(clev[,i])) {
names_num <- c(names_num, i)
}
}
clev_numeric <- clev[,names_num]
clev_cor <- cor(clev_numeric)
which(clev_cor > 0.5 & clev_cor < 1, arr.ind = TRUE) row col
V34 10 2
V37 12 2
V15 5 4
V14 4 5
V31 7 6
V29 6 7
V32 8 7
V31 7 8
V10 2 10
V37 12 11
V10 2 12
V35 11 12which(-clev_cor > 0.5 & -clev_cor < 1, arr.ind = TRUE) row col2.2 Modification of values
par(mfrow = c(2,3))
for(i in 1:length(clev_numeric)){
qqnorm(clev_numeric[,i], main = c("Q-Q Plot: ", names(clev_numeric)[i]))
qqline(clev_numeric[,i], col=2)
}
Boxcox
par(mfrow = c(2,3))
for(i in 1:length(clev_numeric)){
boxcox(lm(clev_numeric[,i]-min(clev_numeric[,i])+1e-6~1),lambda = seq(-1, 1.5, by=0.1), xlab = c(names(clev_numeric))[i])
}
par(mfrow = c(1,3))
#we treat them as special cases
aux <- clev_numeric[,"V14"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))
aux <- clev_numeric[,"V15"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))aux <- clev_numeric[,"V40"]
aux <- aux[aux !=0]
boxcox(lm(aux~1),lambda = seq(-2, 1.5, by=0.1))
clev_sqrt <- c("V10", "V12", "V31", "V43")
clev_sqrt_especial <- c("V14", "V40")
clev_box <- clev_numeric
#box-cox transformation
for (i in 1:ncol(clev_box)){
if (names(clev_box)[i] %in% clev_sqrt) {
clev_box[,i] <- 2*sqrt(clev_box[,i]-min(clev_box[,i])+1e-6)
} else if (names(clev_box)[i] %in% clev_sqrt_especial){
clev_box[,i] <- 2*sqrt(clev_box[,i])
}
}
clev_box <- data.frame(scale(clev_box, scale = T)) #standarized
names(clev_box) <- names(clev_numeric)
clev_box$V14 <- clev_box$V14 - min(clev_box$V14)
clev_box$V15 <- clev_box$V15 - min(clev_box$V15)
clev_box$V40 <- clev_box$V40 - min(clev_box$V40)par(mfrow = c(2,3))
for(i in 1:ncol(clev_box)){
qqnorm(clev_box[,i], main = c("Q-Q Plot: ", names(clev_box)[i]))
qqline(clev_box[,i], col=2)
}
For further models, normalisation will be useful.
And this is the final dataset after processing it.
clev[, names(clev_box)] <- clev_box[,names(clev_box)]
summary(clev) V3 V4 V9 V10 V11
Min. :29.00 0: 91 1: 22 Min. : 94.0 0:108
1st Qu.:48.00 1:191 2: 43 1st Qu.:120.0 1:174
Median :55.00 3: 84 Median :130.0
Mean :54.41 4:133 Mean :131.6
3rd Qu.:61.00 3rd Qu.:140.0
Max. :77.00 Max. :200.0
V12 V14 V15 V16 V18
Min. :126.0 Min. : 0.00 Min. : 0.00 0:240 0:107
1st Qu.:213.0 1st Qu.: 0.00 1st Qu.: 0.00 1: 42 1:175
Median :244.0 Median :10.00 Median :15.00
Mean :249.1 Mean :16.64 Mean :15.08
3rd Qu.:277.0 3rd Qu.:30.00 3rd Qu.:30.00
Max. :564.0 Max. :99.00 Max. :54.00
V19 V20 V21 V22 V23 V24
0:138 1 :38 4 : 13 81:67 0:271 0:186
1: 2 12 :33 15 : 13 82:96 1: 11 1: 95
2:142 2 :31 20 : 13 83:86 2: 1
8 :30 13 : 12 84:33
6 :26 16 : 12
11 :25 21 : 12
(Other):99 (Other):207
V25 V26 V27 V29 V31
0:211 0:252 0:248 Min. : 1.800 Min. : 3.000
1: 71 1: 30 1: 34 1st Qu.: 6.500 1st Qu.: 7.000
Median : 8.500 Median : 9.000
Mean : 8.418 Mean : 9.754
3rd Qu.:10.075 3rd Qu.:12.000
Max. :15.000 Max. :18.000
V32 V33 V34
Min. : 71.0 Min. : 40.00 Min. : 84.0
1st Qu.:133.2 1st Qu.: 65.00 1st Qu.:154.0
Median :153.5 Median : 74.00 Median :168.0
Mean :149.8 Mean : 75.12 Mean :168.1
3rd Qu.:165.8 3rd Qu.: 84.00 3rd Qu.:184.0
Max. :202.0 Max. :119.00 Max. :232.0
V35 V37 V38 V39
Min. : 26.00 Min. : 50.00 0:190 0:276
1st Qu.: 70.00 1st Qu.: 80.00 1: 92 1: 6
Median : 80.00 Median : 85.00
Mean : 78.74 Mean : 84.95
3rd Qu.: 85.00 3rd Qu.: 90.00
Max. :120.00 Max. :110.00
V40 V41 V43 V44 V51
Min. :0.000 1:135 Min. : 24.0 Min. :0.0000 1: 2
1st Qu.:0.000 2:129 1st Qu.: 92.0 1st Qu.:0.0000 3:159
Median :0.800 3: 18 Median :118.0 Median :0.0000 6: 14
Mean :1.027 Mean :123.6 Mean :0.6702 7:107
3rd Qu.:1.600 3rd Qu.:152.8 3rd Qu.:1.0000
Max. :6.200 Max. :270.0 Max. :3.0000
V56 V58 V59 V60 V61 V63 V65
5 : 14 0:157 1:270 1:242 1:224 1:238 1:236
21 : 14 1: 50 2: 12 2: 40 2: 58 2: 44 2: 46
14 : 12 2: 31
1 : 11 3: 32
17 : 11 4: 12
30 : 11
(Other):209
V67 V68
1:233 1:246
2: 49 2: 36
2.3. Feature extraction
Separe train and test data, seed for reproducibility.
set.seed(2000)
n <- nrow(clev)
train.lenght <- round(2*n/3)
clev <- clev[sample(n),]
train <- clev[1:train.lenght,]
test <- clev[(train.lenght+1):n,]names_num <- c()
for(i in 1:ncol(train)){
if (!is.factor(train[,i])) {
names_num <- c(names_num, i)
}
}
train_num <- train[,names_num]Extract PCA features.
pca <- princomp(train_num)
screeplot(pca)
summary(pca)Importance of components:
Comp.1 Comp.2 Comp.3 Comp.4
Standard deviation 1.8444732 1.5577035 1.3875402 1.21213461
Proportion of Variance 0.2297381 0.1638543 0.1300108 0.09921788
Cumulative Proportion 0.2297381 0.3935924 0.5236032 0.62282104
Comp.5 Comp.6 Comp.7 Comp.8
Standard deviation 1.01035405 0.97964763 0.85803949 0.7938423
Proportion of Variance 0.06893431 0.06480791 0.04971676 0.0425556
Cumulative Proportion 0.69175535 0.75656326 0.80628002 0.8488356
Comp.9 Comp.10 Comp.11
Standard deviation 0.73133601 0.69242679 0.62465621
Proportion of Variance 0.03611787 0.03237695 0.02634938
Cumulative Proportion 0.88495350 0.91733045 0.94367983
Comp.12 Comp.13 Comp.14
Standard deviation 0.57654420 0.50756673 0.43609672
Proportion of Variance 0.02244675 0.01739701 0.01284263
Cumulative Proportion 0.96612658 0.98352359 0.99636622
Comp.15
Standard deviation 0.231971913
Proportion of Variance 0.003633784
Cumulative Proportion 1.000000000
Fp <- pca$scores
Gs <- pca$loadings
Fs <- Fp %*% diag(1/pca$sdev)
Gp <- Gs %*% diag(pca$sdev) * 2
col.class <- as.numeric(train$V58)
col.class[col.class==1] <- "red"
col.class[col.class==2] <- "green"
col.class[col.class==3] <- "blue"
col.class[col.class==4] <- "yellow"
col.class[col.class==5] <- "purple"
plot(Fs[,1], Fs[,2], asp=1, col = col.class, xlab = "First principal component", ylab = "Second principal component")
arrows(rep(0,dim(Gs)[1]),rep(0,dim(Gs)[1]), Gp[,1], Gp[,2])text(Gp[,1], Gp[,2], names(train_num), col = "black")
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))
V1, V59, V57 creates many problems
problematic <- c("V57", "V59")
train <- train[, !(names(train) %in% problematic)]
fda <- lda(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)
#plot(fda)
loadings <- predict(fda)$x
plot(loadings, col = col.class)
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))
train$LD1 <- loadings[,1]
train$LD2 <- loadings[,2]
train$LD3 <- loadings[,3]
train$LD4 <- loadings[,4]
fda_test <- predict(fda, newdata = test)
test$LD1 <- fda_test$x[,1]
test$LD2 <- fda_test$x[,2]
test$LD3 <- fda_test$x[,3]
test$LD4 <- fda_test$x[,4]Análisis por correspondencias
ac <- mjca(clev[,names(clev) %in% factores], lambda="Burt")
plot(ac, main="MCA biplot of Burt matrix with data")
ac_ind <- mjca(train[,names(train) %in% factores], lambda="indicator", reti = T)
plot(ac_ind$rowcoord, col = col.class)
legend("bottomright", fill=c("red","green", "blue", "yellow", "purple"), legend=c('0','1','2','3', '4'))
mca.features <- ac_ind$rowcoord3. Resampling protocol
Principal utility function for doing cross-validation.
First we separate in train and test.
Create the CV function for any model with associated predict and update functions.
cross.validation <- function(data, target, model, times, nfolds, need.class){
set.seed(0202)
CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)
err.total <- c()
for (i in 1:times){
err.onetime <- c()
for (j in 1:nfolds){
print(paste0("Fold: ", j))
val <- unlist(CV.folds[[i]][[j]])
tr <- data[-val,]
va <- data[val,]
model <- update(model, data=tr)
pred <- predict(model, newdata=va)
if (need.class){
pred <- pred$class
}
err.table <- table(True=target[val], Predicted=pred)
err <- 1-sum(diag(err.table))/sum(err.table)
err.onetime <- c(err.onetime, err)
}
err.total <- c(err.total, mean(err.onetime))
print(paste0("Iteration ", i, ", mean error: ", mean(err.onetime)))
}
mean(err.total)
}Cross-validation for k-Nearest Neighbour
cross.validation.knn <- function(data, target, times, nfolds, K){
set.seed(0202)
CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)
err.total <- c()
for (i in 1:times){
err.onetime <- c()
for (j in 1:nfolds){
val <- unlist(CV.folds[[i]][[j]])
tr <- data[-val,]
va <- data[val,]
pred <- knn(tr, va, target[-val], k = K)
err.table <- table(True=target[val], Predicted=pred)
err <- 1-sum(diag(err.table))/sum(err.table)
err.onetime <- c(err.onetime, err)
}
err.total <- c(err.total, mean(err.onetime))
}
mean(err.total)
}Cross-validation Naive-Bayes
cross.validation.naive <- function(data, target, model, times, nfolds){
set.seed(0202)
CV.folds <- generateCVRuns(target, ntimes=times, nfold=nfolds, stratified=TRUE)
err.total <- c()
for (i in 1:times){
err.onetime <- c()
for (j in 1:nfolds){
print(paste0("Fold: ", j))
val <- unlist(CV.folds[[i]][[j]])
tr <- data[-val,]
va <- data[val,]
model <- naiveBayes(V58~.-LD1-LD2-LD3-LD4, data=tr)
pred <- predict(model, newdata=va)
err.table <- table(True=target[val], Predicted=pred)
err <- 1-sum(diag(err.table))/sum(err.table)
err.onetime <- c(err.onetime, err)
}
err.total <- c(err.total, mean(err.onetime))
print(paste0("Iteration ", i, ", mean error: ", mean(err.onetime)))
}
mean(err.total)
}4. Models
The models we are going to use are: - LDA - QDA - RDA - k-NN - Naïve Bayes - GLM
rda.model <- rda(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V20+V21+V22+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)
naive.model <- naiveBayes(V58~V3+V4+V9+V10+V11+V12+V14+V15+V16+V18+V19+V20+V21+V22+V23+V24+V25+V26+V27+V29+V31+V32+V33+V34+V35+V37+V38+V39+V40+V41+V43+V44+V51+V56+V60+V61+V63+V65+V67+V68, data=train)cross.validation(train, train$V58, rda.model, 10, 10, T)rda.model.fda <- rda(V58~.,data=train)cross.validation(train, train$V58, rda.model.fda, 10, 10, T)cross.validation.naive(train, train$V58, naive.model, 10, 10)err <- c()
for (k in 1:20){
err <- c(err, cross.validation.knn(train, train$V58, 10,10, k))
}plot(err, type = "l")
err [1] 0.1519883 0.1965205 0.1950585 0.2004678 0.2073977 0.2152924
[7] 0.2184795 0.2276023 0.2313450 0.2323392 0.2419883 0.2456433
[13] 0.2642398 0.2690351 0.2743275 0.2791228 0.2871053 0.2955556
[19] 0.3003801 0.3050877cross.validation.knn(train, train$V58, 10, 10, 1)
multinomial.model <- multinom(V58~., data=train)
cross.validation(train, train$V58, multinomial.model, 10, 10, F)
multinomial.model.step <- step(multinomial.model)
cross.validation(train, train$V58, multinomial.model.step, 10, 10, F)
multinomial.model.noFDA <- multinom(V58~.-LD1-LD2-LD3-LD4, data=train)
cross.validation(train, train$V58, multinomial.model.noFDA, 10, 10, F)
multinomial.model.noFDA.step <- step(multinomial.model.noFDA)
cross.validation(train, train$V58, multinomial.model.noFDA.step, 10, 10, F)Test error
rda.model <- update(rda.model.fda, data=train)
pred.test <- predict(rda.model.fda, test)
pred.test <- pred.test$class
(err.table <- table(True=test$V58, Pred=pred.test)) Pred
True 0 1 2 3 4
0 40 5 0 0 0
1 6 13 1 0 2
2 0 2 6 4 1
3 0 0 1 9 0
4 0 0 0 1 3(err.test <- 1-sum(diag(err.table))/sum(err.table))[1] 0.2446809Parte 2
We will now use model on a dataset not yet investigated: Hungary.
hung <- read.table("../data/processed.hungarian.data", header=F, sep=',', na.strings="?")
summary(hung) V1 V2 V3 V4
Min. :28.00 Min. :0.0000 Min. :1.000 Min. : 92.0
1st Qu.:42.00 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.:120.0
Median :49.00 Median :1.0000 Median :3.000 Median :130.0
Mean :47.83 Mean :0.7245 Mean :2.983 Mean :132.6
3rd Qu.:54.00 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:140.0
Max. :66.00 Max. :1.0000 Max. :4.000 Max. :200.0
NA's :1
V5 V6 V7
Min. : 85.0 Min. :0.00000 Min. :0.0000
1st Qu.:209.0 1st Qu.:0.00000 1st Qu.:0.0000
Median :243.0 Median :0.00000 Median :0.0000
Mean :250.8 Mean :0.06993 Mean :0.2184
3rd Qu.:282.5 3rd Qu.:0.00000 3rd Qu.:0.0000
Max. :603.0 Max. :1.00000 Max. :2.0000
NA's :23 NA's :8 NA's :1
V8 V9 V10
Min. : 82.0 Min. :0.0000 Min. :0.0000
1st Qu.:122.0 1st Qu.:0.0000 1st Qu.:0.0000
Median :140.0 Median :0.0000 Median :0.0000
Mean :139.1 Mean :0.3038 Mean :0.5861
3rd Qu.:155.0 3rd Qu.:1.0000 3rd Qu.:1.0000
Max. :190.0 Max. :1.0000 Max. :5.0000
NA's :1 NA's :1
V11 V12 V13 V14
Min. :1.000 Min. :0 Min. :3.000 Min. :0.0000
1st Qu.:2.000 1st Qu.:0 1st Qu.:5.250 1st Qu.:0.0000
Median :2.000 Median :0 Median :6.000 Median :0.0000
Mean :1.894 Mean :0 Mean :5.643 Mean :0.3605
3rd Qu.:2.000 3rd Qu.:0 3rd Qu.:7.000 3rd Qu.:1.0000
Max. :3.000 Max. :0 Max. :7.000 Max. :1.0000
NA's :190 NA's :291 NA's :266 Missing values treatment
rm.var <- c("V11", "V12", "V13")
hung <- hung[,!(names(hung) %in% rm.var)]hung[is.na(hung)] <- -9
var.imputable <- c("V4", "V5", "V6", "V7", "V8", "V9")
for (nom in var.imputable){
hung[, nom][hung[, nom] == -9] <- NA
variable <- hung[, nom]
hung[, nom] <- knn.imputation(hung, variable, nom, 7)
}var.factor <- c("V2", "V3", "V6", "V7", "V9", "V14")
for (nom in var.factor){
hung[, nom] <- as.factor(as.character(hung[,nom]))
}summary(hung)
V1 V2 V3 V4 V5
Min. :28.00 0: 81 1: 11 Min. : 27.0 Min. : 4.0
1st Qu.:42.00 1:213 2:106 1st Qu.:120.0 1st Qu.:198.0
Median :49.00 3: 54 Median :130.0 Median :237.0
Mean :47.83 4:123 Mean :132.2 Mean :236.6
3rd Qu.:54.00 3rd Qu.:140.0 3rd Qu.:277.0
Max. :66.00 Max. :200.0 Max. :603.0
V6 V7 V8 V9 V10 V14
0:266 0:235 Min. : 54.0 0:204 Min. :0.0000 0:188
1: 28 1: 53 1st Qu.:122.0 1: 89 1st Qu.:0.0000 1:106
2: 6 Median :140.0 2: 1 Median :0.0000
Mean :138.8 Mean :0.5861
3rd Qu.:155.0 3rd Qu.:1.0000
Max. :190.0 Max. :5.0000 Gaussianity treatment
names_num <- c()
for(i in 1:ncol(hung)){
if (!is.factor(hung[,i])) {
names_num <- c(names_num, i)
}
}
hung_numeric <- hung[,names_num]
hung_cor <- cor(hung_numeric)
which(hung_cor > 0.5 & hung_cor < 1, arr.ind = TRUE) row colwhich(-hung_cor > 0.5 & -hung_cor < 1, arr.ind = TRUE) row colpar(mfrow = c(2,3))
for(i in 1:length(hung_numeric)){
qqnorm(hung_numeric[,i], main = c("Q-Q Plot: ", names(hung_numeric)[i]))
qqline(hung_numeric[,i], col=2)
}
Boxcox
par(mfrow = c(2,3))
for(i in 1:length(hung_numeric)){
boxcox(lm(hung_numeric[,i]-min(hung_numeric[,i])+1e-6~1),lambda = seq(-1, 1.5, by=0.1), xlab = c(names(hung_numeric))[i])
}
We see var 10 requires a logaritmic transformation.
qqnorm(hung[,"V10"])
qqline(hung[,"V10"], col=2)
Feature extraction
Separe train and test data, seed for reproducibility.
set.seed(2000)
n <- nrow(hung)
train.lenght <- round(2*n/3)
hung <- hung[sample(n),]
train <- hung[1:train.lenght,]
test <- hung[(train.lenght+1):n,]names_num <- c()
for(i in 1:ncol(train)){
if (!is.factor(train[,i])) {
names_num <- c(names_num, i)
}
}
train_num <- train[,names_num]Extract PCA features.
pca <- princomp(train_num)
screeplot(pca)
summary(pca)Importance of components:
Comp.1 Comp.2 Comp.3
Standard deviation 86.0161788 24.47912423 17.80405417
Proportion of Variance 0.8849996 0.07167612 0.03791583
Cumulative Proportion 0.8849996 0.95667568 0.99459151
Comp.4 Comp.5
Standard deviation 6.709559731 0.4448637716
Proportion of Variance 0.005384815 0.0000236721
Cumulative Proportion 0.999976328 1.0000000000
Fp <- pca$scores
Gs <- pca$loadings
Fs <- Fp %*% diag(1/pca$sdev)
Gp <- Gs %*% diag(pca$sdev) *0.3
col.class <- as.numeric(train$V14)
col.class[col.class==0] <- "red"
col.class[col.class==1] <- "blue"
plot(Fs[,1], Fs[,2], asp=1, col = col.class, xlab = "First principal component", ylab = "Second principal component")
arrows(rep(0,dim(Gs)[1]),rep(0,dim(Gs)[1]), Gp[,1], Gp[,2])text(Gp[,1], Gp[,2], names(train_num), col = "black")
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))
fda <- lda(V14~., data=train)
#plot(fda)
loadings <- predict(fda)$x
plot(loadings, col = col.class)
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))
train$LD1 <- loadings[,1]
fda_test <- predict(fda, newdata = test)
test$LD1 <- fda_test$x[,1]Análisis por correspondencias
ac <- mjca(hung[,names(hung) %in% var.factor], lambda="Burt")
plot(ac, main="MCA biplot of Burt matrix with data")
ac_ind <- mjca(train[,names(train) %in% var.factor], lambda="indicator", reti = T)
plot(ac_ind$rowcoord, col = col.class)
legend("bottomright", fill=c("red","blue"), legend=c('0','1'))
mca.features <- ac_ind$rowcoordLS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCiMgTWFjaGluZSBsZWFybmluZyBwcm9qZWN0CiMjIEF1dGhvcnM6IEpvc2UgUMOpcmV6IENhbm8gJiDDgWx2YXJvIFJpYm90IEJhcnJhZG8KCiMjIyAwLiBMaWJyYXJpZXMKCmBgYHtyfQppbnN0YWxsLnBhY2thZ2VzKCJrbGFSIikKaW5zdGFsbC5wYWNrYWdlcygiVHVuZVBhcmV0byIpCmluc3RhbGwucGFja2FnZXMoInJnbCIpCmluc3RhbGwucGFja2FnZXMoImdsbW5ldCIpCmluc3RhbGwucGFja2FnZXMoImNhIikKYGBgCgpgYGB7cn0KIyBMREEvIFFEQQpsaWJyYXJ5KE1BU1MpCgojIFJEQQpsaWJyYXJ5KGtsYVIpCgojIE11bHRpbm9taWFsCmxpYnJhcnkobm5ldCkKCiMgQ3Jvc3MtVmFsaWRhdGlvbgpsaWJyYXJ5KFR1bmVQYXJldG8pCgojIE5haXZlIEJheWVzCmxpYnJhcnkoZTEwNzEpCgojIGstTk4KbGlicmFyeShjbGFzcykKCiMgM2QKbGlicmFyeShyZ2wpCgojIExBU1NPCmxpYnJhcnkoTWF0cml4KQpsaWJyYXJ5KGdsbW5ldCkKCiMgQ29ycmVzcG9uZGVuY2UgYW5hbHlzaXMKbGlicmFyeShjYSkKYGBgCgoKIyMjIDEuIFJlYWQgZGF0YQoKYGBge3J9CnNldC5zZWVkKDIxMDUpCnNldHdkKCIuLi9kYXRhIikKY2xldiA8LSByZWFkLmNzdigiY2xldmVsYW5kLmNzdiIsIGhlYWRlcj1GKQpodW5nIDwtIHJlYWQuY3N2KCJodW5nYXJpYW4uY3N2IiwgaGVhZGVyPUYpCnZhIDwtIHJlYWQuY3N2KCJsb25nLWJlYWNoLXZhLmNzdiIsIGhlYWRlcj1GKQpzd2l0eiA8LSByZWFkLmNzdigic3dpdHplcmxhbmQuY3N2IiwgaGVhZGVyPUYpCgpjbGV2JGxvY2F0aW9uIDwtICJjbGV2ZWxhbmQiCmh1bmckbG9jYXRpb24gPC0gImh1bmdhcmlhbiIKdmEkbG9jYXRpb24gPC0gImxvbmctYmVhY2gtdmEiCnN3aXR6JGxvY2F0aW9uIDwtICJzd2l0emVybGFuZCIKCmhlYXJ0MSA8LSByYmluZChjbGV2LCBodW5nKQpoZWFydDIgPC0gcmJpbmQodmEsIHN3aXR6KQpoZWFydCA8LSByYmluZChoZWFydDEsIGhlYXJ0MikKaGVhZChoZWFydCkKYGBgCgoKIyMjIDIuIFByZXByb2Nlc3MgZGF0YQpXZSBhcHBseSBjbHVzdGVyaW5nIGFuZCBzZXZlcmFsIHBsb3R0aW5nIHRlY2huaXF1ZXMgdG8gaGF2ZSBhbiBpZGVhIG9mIHRoZSBkYXRhc2V0LiBJbiBjYXNlIHRoZXJlIGFyZSBOQXMgd2Ugd2lsbCB1c2Ugay1OTiBmb3IgaW1wdXRhdGlvbi4gCgpUbyBleHRyYWN0IG5ldyBmZWF0dXJlcyB3ZSB3aWxsIGFwcGx5IFBDQSBhbmQgRkRBIGFuZCBrZWVwIHRoaXMgbmV3IGZlYXR1cmVzIGFuZCBjb21wb25lbnRzIGFwYXJ0LgoKYGBge3J9CiMgSXQgc2F5cyB3aGljaCBjb2x1bW5zIGFyZSBhbGwgbWlzc2luZ3MKIyBUaGUgaW5kZXggYXJlIHJldHVybmVkIGluIG5lZ2F0aXZlIHRvIGVsaW1pbmF0ZSB0aGVtCm5hLmNvbHVtbnMgPC0gZnVuY3Rpb24oZGQpewogIHJtbGlzdCA8LSBjKCkKICBmb3IgKGkgaW4gMTpuY29sKGRkKSl7CiAgICBpZiAobWluKGRkWyxpXSkgPT0gbWF4KGRkWyxpXSkgJiBtaW4oZGRbLGldKT09LTkpewogICAgICBybWxpc3QgPC0gYyhybWxpc3QsIGkpCiAgICB9CiAgfQogIC1ybWxpc3QKfQoKY2xldiA8LSBjbGV2WyxuYS5jb2x1bW5zKGNsZXYpXQoKIyBSZXR1cm5zIGNvbHVtbnMgd2l0aCBtb3JlIE5BIHRoYW4gYSBnaXZlbiB0aHJlc2hvbGQsIGFsc28gaW4gbmVnYXRpdmUKbXVjaC5uYS5jb2xzIDwtIGZ1bmN0aW9uKGRkLCB0aHJlc2hvbGQpewogIHJtbGlzdCA8LSBjKCkKICBmb3IgKGkgaW4gMTpuY29sKGRkKSl7CiAgICBpZiAoc3VtKGRkWyxpXT09LTkpID4gdGhyZXNob2xkKXsKICAgICAgcm1saXN0IDwtIGMocm1saXN0LCBpKQogICAgfQogIH0gICAKICAtcm1saXN0ICAKfQoKY2xldiA8LSBjbGV2WywgbXVjaC5uYS5jb2xzKGNsZXYsIDYwKV0KCiMgQXBwbGllcyBrLW5lYXJlc3QgbmVpZ2hib3VyIGltcHV0YXRpb24gZm9yIGEgZ2l2ZW4gdmFyaWFibGUKa25uLmltcHV0YXRpb24gPSBmdW5jdGlvbiAoZGQsIHZhcmlhYmxlLCB2YXJuYW1lLCBrKQp7ICAKICBhdXggPSBzdWJzZXQgKGRkLCBzZWxlY3QgPSBuYW1lcyhkZClbbmFtZXMoZGQpICE9IHZhcm5hbWVdKQogIGF1eDEgPSBhdXhbIWlzLm5hKHZhcmlhYmxlKSxdCiAgYXV4MiA9IGF1eFtpcy5uYSh2YXJpYWJsZSksXQoKICAjIE5laXRoZXIgb2YgYXV4MSwgYXV4MiBjYW4gY29udGFpbiBOQXMKICBrbm4uaW5jID0ga25uIChhdXgxLGF1eDIsIHZhcmlhYmxlWyFpcy5uYSh2YXJpYWJsZSldLCBrKQogIHZhcmlhYmxlW2lzLm5hKHZhcmlhYmxlKV0gPSBrbm4uaW5jCiAgdmFyaWFibGUKfQoKIyBUaGlzIGFyZSB0aGUgdmFyaWFibGVzIHdoaWNoIHZhbHVlcyB3aGVyZSBzdWJzdGl0dXRlZCBieSBkdW1teSB2YWx1ZXMuCmR1bW15IDwtIGMoIlYxIiwgIlYyIiwgIlYzNiIsICJWNjkiLCAiVjcwIiwgIlY3MSIsICJWNzIiLCAiVjczIiwgIlYyOCIsICJsb2NhdGlvbiIpCmNsZXYgPC0gY2xldlssIShuYW1lcyhjbGV2KSAlaW4lIGR1bW15KV0KCiMga25uIGltcHV0YXRpb24gZm9yIGNsZXYKbmEubmFtZXMgPC0gbmFtZXMoY2xldilbLW11Y2gubmEuY29scyhjbGV2LCAwKV0KZm9yIChuYW1lIGluIG5hLm5hbWVzKXsKICBjbGV2WywgbmFtZV1bY2xldlssIG5hbWVdID09IC05XSA8LSBOQQogIGNsZXZbLCBuYW1lXSA8LSBrbm4uaW1wdXRhdGlvbihjbGV2LCBjbGV2WyxuYW1lXSwgbmFtZSwgNykKfQpgYGAKCk5vdywgd2UgYW5hbHlzZSB0aGUgY29ycmVsYXRpb25zIGFtb25nIHZhcmlhYmxlcyBzaW5jZSBpdCB3aWxsIGhhdmUgaW1wYWN0IG9uIGxhdGVyIG1vZGVscy4KCmBgYHtyfQpjb3JyLmZhY3RvcnMgPC0gY29yKGNsZXYpCndoaWNoKGFicyhjb3JyLmZhY3RvcnMpLWRpYWcoZGlhZyhjb3JyLmZhY3RvcnMpKT4wLjksIGFyci5pbmQ9VCkKCnJtLmNvcnJlbGF0ZWQgPC0gYygiVjU3IiwgIlY1NSIpCmNsZXYgPC0gY2xldlssIShuYW1lcyhjbGV2KSAlaW4lIHJtLmNvcnJlbGF0ZWQpXQoKZmFjdG9yZXMgPC0gYygiVjU4IiwgIlY0IiwgIlY5IiwgIlYxNiIsICJWMTgiLCAiVjE5IiwgIlYyMCIsICJWMjEiLCAiVjIyIiwgIlYyMyIsICJWMjQiLCAiVjI1IiwgIlYyNiIsICJWMjciLCAiVjM4IiwgIlYzOSIsICJWNDEiLCAiVjUxIiwgIlY1NiIsICJWMTEiLCAiVjU5IiwgIlY2MCIsICJWNjEiLCAiVjYzIiwgIlY2NSIsICJWNjciLCAiVjY4IikKZm9yIChmIGluIGZhY3RvcmVzKXsKICBjbGV2WyxmXSA8LSBhcy5mYWN0b3IoY2xldlssZl0pCn0KCiMgRHVtbXkgbGV2ZWwgZ2V0cyByZXBsYWNlZApjbGV2JFYyNVtjbGV2JFYyNSA9PSAyXSA8LSAxCmNsZXYkVjI1IDwtIGRyb3BsZXZlbHMoY2xldiRWMjUpCmxldmVscyhjbGV2JFYyNSkKYGBgCgpUaGlzIGlzIHRoZSBkYXRhc2V0IGFmdGVyIHRoZSB0cmVhdG1lbnQgZm9yIG1pc3NpbmdzLgoKYGBge3J9CnN1bW1hcnkoY2xldikKYGBgCgoKIyMjIyAyLjEgVmlzdWFsaXphdGlvbnMKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKCFpcy5mYWN0b3IoY2xldlssaV0pKSB7CiAgICBoaXN0KGNsZXZbLGldLCBtYWluID0gbmFtZXMoY2xldilbaV0sIHhsYWI9IlZhbHVlcyIpCiAgfQp9CmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKCFpcy5mYWN0b3IoY2xldlssaV0pKSB7CiAgICBib3hwbG90KGNsZXZbLGldLCB4bGFiID0gbmFtZXMoY2xldilbaV0pCiAgfQp9CmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDMsMykpCmZvcihpIGluIDE6bmNvbChjbGV2KSl7CiAgaWYgKGlzLmZhY3RvcihjbGV2WyxpXSkpIHsKICAgIGhpc3QoYXMubnVtZXJpYyhhcy5jaGFyYWN0ZXIoY2xldlssaV0pKSwgbWFpbiA9IG5hbWVzKGNsZXYpW2ldLCB4bGFiPSJWYWx1ZXMiKQogIH0KfQpgYGAKCgpgYGB7cn0KbmFtZXNfbnVtIDwtIGMoKQpmb3IoaSBpbiAxOm5jb2woY2xldikpewogIGlmICghaXMuZmFjdG9yKGNsZXZbLGldKSkgewogICAgbmFtZXNfbnVtIDwtIGMobmFtZXNfbnVtLCBpKQogIH0KfQpjbGV2X251bWVyaWMgPC0gY2xldlssbmFtZXNfbnVtXQoKY2xldl9jb3IgPC0gY29yKGNsZXZfbnVtZXJpYykKd2hpY2goY2xldl9jb3IgPiAwLjUgJiBjbGV2X2NvciA8IDEsIGFyci5pbmQgPSBUUlVFKQp3aGljaCgtY2xldl9jb3IgPiAwLjUgJiAtY2xldl9jb3IgPCAxLCBhcnIuaW5kID0gVFJVRSkKYGBgCgoKIyMjIyAyLjIgTW9kaWZpY2F0aW9uIG9mIHZhbHVlcwoKYGBge3J9CnBhcihtZnJvdyA9IGMoMiwzKSkKZm9yKGkgaW4gMTpsZW5ndGgoY2xldl9udW1lcmljKSl7CiAgICBxcW5vcm0oY2xldl9udW1lcmljWyxpXSwgbWFpbiA9IGMoIlEtUSBQbG90OiAiLCBuYW1lcyhjbGV2X251bWVyaWMpW2ldKSkKICAgIHFxbGluZShjbGV2X251bWVyaWNbLGldLCBjb2w9MikKfQpgYGAKCgpCb3hjb3gKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bGVuZ3RoKGNsZXZfbnVtZXJpYykpewogIGJveGNveChsbShjbGV2X251bWVyaWNbLGldLW1pbihjbGV2X251bWVyaWNbLGldKSsxZS02fjEpLGxhbWJkYSA9IHNlcSgtMSwgMS41LCBieT0wLjEpLCB4bGFiID0gYyhuYW1lcyhjbGV2X251bWVyaWMpKVtpXSkKfQpgYGAKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDEsMykpCiN3ZSB0cmVhdCB0aGVtIGFzIHNwZWNpYWwgY2FzZXMKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlYxNCJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlYxNSJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYXV4IDwtIGNsZXZfbnVtZXJpY1ssIlY0MCJdCmF1eCA8LSBhdXhbYXV4ICE9MF0KYm94Y294KGxtKGF1eH4xKSxsYW1iZGEgPSBzZXEoLTIsIDEuNSwgYnk9MC4xKSkKYGBgCgoKYGBge3J9CmNsZXZfc3FydCA8LSBjKCJWMTAiLCAiVjEyIiwgIlYzMSIsICJWNDMiKQpjbGV2X3NxcnRfZXNwZWNpYWwgPC0gYygiVjE0IiwgIlY0MCIpCmNsZXZfYm94IDwtIGNsZXZfbnVtZXJpYwojYm94LWNveCB0cmFuc2Zvcm1hdGlvbgpmb3IgKGkgaW4gMTpuY29sKGNsZXZfYm94KSl7CiAgaWYgKG5hbWVzKGNsZXZfYm94KVtpXSAlaW4lIGNsZXZfc3FydCkgewogICAgY2xldl9ib3hbLGldIDwtIDIqc3FydChjbGV2X2JveFssaV0tbWluKGNsZXZfYm94WyxpXSkrMWUtNikKICB9IGVsc2UgaWYgKG5hbWVzKGNsZXZfYm94KVtpXSAlaW4lIGNsZXZfc3FydF9lc3BlY2lhbCl7CiAgICBjbGV2X2JveFssaV0gPC0gMipzcXJ0KGNsZXZfYm94WyxpXSkKICB9Cn0KY2xldl9ib3ggPC0gZGF0YS5mcmFtZShzY2FsZShjbGV2X2JveCwgc2NhbGUgPSBUKSkgI3N0YW5kYXJpemVkCm5hbWVzKGNsZXZfYm94KSA8LSBuYW1lcyhjbGV2X251bWVyaWMpCmNsZXZfYm94JFYxNCA8LSBjbGV2X2JveCRWMTQgLSBtaW4oY2xldl9ib3gkVjE0KQpjbGV2X2JveCRWMTUgPC0gY2xldl9ib3gkVjE1IC0gbWluKGNsZXZfYm94JFYxNSkKY2xldl9ib3gkVjQwIDwtIGNsZXZfYm94JFY0MCAtIG1pbihjbGV2X2JveCRWNDApCmBgYAoKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bmNvbChjbGV2X2JveCkpewogICAgcXFub3JtKGNsZXZfYm94WyxpXSwgbWFpbiA9IGMoIlEtUSBQbG90OiAiLCBuYW1lcyhjbGV2X2JveClbaV0pKQogICAgcXFsaW5lKGNsZXZfYm94WyxpXSwgY29sPTIpCn0KYGBgCgpGb3IgZnVydGhlciBtb2RlbHMsIG5vcm1hbGlzYXRpb24gd2lsbCBiZSB1c2VmdWwuCgpgYGB7cn0KY2xldl9ib3ggPC0gYXBwbHkoY2xldl9ib3gsIDIsIHNjYWxlKQpgYGAKCgpBbmQgdGhpcyBpcyB0aGUgZmluYWwgZGF0YXNldCBhZnRlciBwcm9jZXNzaW5nIGl0LgoKYGBge3J9CmNsZXZbLCBuYW1lcyhjbGV2X2JveCldIDwtIGNsZXZfYm94WyxuYW1lcyhjbGV2X2JveCldCnN1bW1hcnkoY2xldikKYGBgCgoKIyMjIyAyLjMuIEZlYXR1cmUgZXh0cmFjdGlvbgoKU2VwYXJlIHRyYWluIGFuZCB0ZXN0IGRhdGEsIHNlZWQgZm9yIHJlcHJvZHVjaWJpbGl0eS4KCmBgYHtyfQpzZXQuc2VlZCgyMDAwKQpuIDwtIG5yb3coY2xldikKdHJhaW4ubGVuZ2h0IDwtIHJvdW5kKDIqbi8zKQoKY2xldiA8LSBjbGV2W3NhbXBsZShuKSxdCnRyYWluIDwtIGNsZXZbMTp0cmFpbi5sZW5naHQsXQp0ZXN0IDwtIGNsZXZbKHRyYWluLmxlbmdodCsxKTpuLF0KYGBgCgpgYGB7cn0KbmFtZXNfbnVtIDwtIGMoKQpmb3IoaSBpbiAxOm5jb2wodHJhaW4pKXsKICBpZiAoIWlzLmZhY3Rvcih0cmFpblssaV0pKSB7CiAgICBuYW1lc19udW0gPC0gYyhuYW1lc19udW0sIGkpCiAgfQp9CnRyYWluX251bSA8LSB0cmFpblssbmFtZXNfbnVtXQpgYGAKCkV4dHJhY3QgUENBIGZlYXR1cmVzLgoKYGBge3J9CnBjYSA8LSBwcmluY29tcCh0cmFpbl9udW0pCnNjcmVlcGxvdChwY2EpCnN1bW1hcnkocGNhKQpgYGAKCgpgYGB7cn0KYmlwbG90KHBjYSkKYGBgCgpgYGB7cn0KRnAgPC0gcGNhJHNjb3JlcwpHcyA8LSBwY2EkbG9hZGluZ3MKCkZzIDwtIEZwICUqJSBkaWFnKDEvcGNhJHNkZXYpCkdwIDwtIEdzICUqJSBkaWFnKHBjYSRzZGV2KSAqIDIKCmNvbC5jbGFzcyA8LSBhcy5udW1lcmljKHRyYWluJFY1OCkKY29sLmNsYXNzW2NvbC5jbGFzcz09MV0gPC0gInJlZCIKY29sLmNsYXNzW2NvbC5jbGFzcz09Ml0gPC0gImdyZWVuIgpjb2wuY2xhc3NbY29sLmNsYXNzPT0zXSA8LSAiYmx1ZSIKY29sLmNsYXNzW2NvbC5jbGFzcz09NF0gPC0gInllbGxvdyIKY29sLmNsYXNzW2NvbC5jbGFzcz09NV0gPC0gInB1cnBsZSIKCnBsb3QoRnNbLDFdLCBGc1ssMl0sIGFzcD0xLCBjb2wgPSBjb2wuY2xhc3MsIHhsYWIgPSAiRmlyc3QgcHJpbmNpcGFsIGNvbXBvbmVudCIsIHlsYWIgPSAiU2Vjb25kIHByaW5jaXBhbCBjb21wb25lbnQiKQphcnJvd3MocmVwKDAsZGltKEdzKVsxXSkscmVwKDAsZGltKEdzKVsxXSksIEdwWywxXSwgR3BbLDJdKQp0ZXh0KEdwWywxXSwgR3BbLDJdLCBuYW1lcyh0cmFpbl9udW0pLCBjb2wgPSAiYmxhY2siKQpsZWdlbmQoImJvdHRvbXJpZ2h0IiwgZmlsbD1jKCJyZWQiLCJncmVlbiIsICJibHVlIiwgInllbGxvdyIsICJwdXJwbGUiKSwgbGVnZW5kPWMoJzAnLCcxJywnMicsJzMnLCAnNCcpKQpgYGAKCgpWMSwgVjU5LCBWNTcgY3JlYXRlcyBtYW55IHByb2JsZW1zCgpgYGB7cn0KcHJvYmxlbWF0aWMgPC0gYygiVjU3IiwgIlY1OSIpCnRyYWluIDwtIHRyYWluWywgIShuYW1lcyh0cmFpbikgJWluJSBwcm9ibGVtYXRpYyldCmZkYSA8LSBsZGEoVjU4flYzK1Y0K1Y5K1YxMCtWMTErVjEyK1YxNCtWMTUrVjE2K1YxOCtWMTkrVjIzK1YyNCtWMjUrVjI2K1YyNytWMjkrVjMxK1YzMitWMzMrVjM0K1YzNStWMzcrVjM4K1YzOStWNDArVjQxK1Y0MytWNDQrVjUxK1Y1NitWNjArVjYxK1Y2MytWNjUrVjY3K1Y2OCwgZGF0YT10cmFpbikKI3Bsb3QoZmRhKQpsb2FkaW5ncyA8LSBwcmVkaWN0KGZkYSkkeApwbG90KGxvYWRpbmdzLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImdyZWVuIiwgImJsdWUiLCAieWVsbG93IiwgInB1cnBsZSIpLCBsZWdlbmQ9YygnMCcsJzEnLCcyJywnMycsICc0JykpCmBgYAoKCmBgYHtyfQp0cmFpbiRMRDEgPC0gbG9hZGluZ3NbLDFdCnRyYWluJExEMiA8LSBsb2FkaW5nc1ssMl0KdHJhaW4kTEQzIDwtIGxvYWRpbmdzWywzXQp0cmFpbiRMRDQgPC0gbG9hZGluZ3NbLDRdCgpmZGFfdGVzdCA8LSBwcmVkaWN0KGZkYSwgbmV3ZGF0YSA9IHRlc3QpCnRlc3QkTEQxIDwtIGZkYV90ZXN0JHhbLDFdCnRlc3QkTEQyIDwtIGZkYV90ZXN0JHhbLDJdCnRlc3QkTEQzIDwtIGZkYV90ZXN0JHhbLDNdCnRlc3QkTEQ0IDwtIGZkYV90ZXN0JHhbLDRdCmBgYAoKCkFuw6FsaXNpcyBwb3IgY29ycmVzcG9uZGVuY2lhcwoKYGBge3J9CmFjIDwtIG1qY2EoY2xldlssbmFtZXMoY2xldikgJWluJSBmYWN0b3Jlc10sIGxhbWJkYT0iQnVydCIpCnBsb3QoYWMsIG1haW49Ik1DQSBiaXBsb3Qgb2YgQnVydCBtYXRyaXggd2l0aCBkYXRhIikKCmFjX2luZCA8LSBtamNhKHRyYWluWyxuYW1lcyh0cmFpbikgJWluJSBmYWN0b3Jlc10sIGxhbWJkYT0iaW5kaWNhdG9yIiwgcmV0aSA9IFQpCnBsb3QoYWNfaW5kJHJvd2Nvb3JkLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImdyZWVuIiwgImJsdWUiLCAieWVsbG93IiwgInB1cnBsZSIpLCBsZWdlbmQ9YygnMCcsJzEnLCcyJywnMycsICc0JykpCgptY2EuZmVhdHVyZXMgPC0gYWNfaW5kJHJvd2Nvb3JkCmBgYAoKCiMjIyAzLiBSZXNhbXBsaW5nIHByb3RvY29sClByaW5jaXBhbCB1dGlsaXR5IGZ1bmN0aW9uIGZvciBkb2luZyBjcm9zcy12YWxpZGF0aW9uLiAKCkZpcnN0IHdlIHNlcGFyYXRlIGluIHRyYWluIGFuZCB0ZXN0LgoKQ3JlYXRlIHRoZSBDViBmdW5jdGlvbiBmb3IgYW55IG1vZGVsIHdpdGggYXNzb2NpYXRlZCBwcmVkaWN0IGFuZCB1cGRhdGUgZnVuY3Rpb25zLgoKYGBge3J9CmNyb3NzLnZhbGlkYXRpb24gPC0gZnVuY3Rpb24oZGF0YSwgdGFyZ2V0LCBtb2RlbCwgdGltZXMsIG5mb2xkcywgbmVlZC5jbGFzcyl7CiAgc2V0LnNlZWQoMDIwMikKICBDVi5mb2xkcyA8LSBnZW5lcmF0ZUNWUnVucyh0YXJnZXQsIG50aW1lcz10aW1lcywgbmZvbGQ9bmZvbGRzLCBzdHJhdGlmaWVkPVRSVUUpCgogIGVyci50b3RhbCA8LSBjKCkKICBmb3IgKGkgaW4gMTp0aW1lcyl7CiAgICBlcnIub25ldGltZSA8LSBjKCkKICAgIGZvciAoaiBpbiAxOm5mb2xkcyl7CiAgICAgIHByaW50KHBhc3RlMCgiRm9sZDogIiwgaikpCiAgICAgIHZhbCA8LSB1bmxpc3QoQ1YuZm9sZHNbW2ldXVtbal1dKQogICAgICAKICAgICAgdHIgPC0gZGF0YVstdmFsLF0KICAgICAgdmEgPC0gZGF0YVt2YWwsXQoKICAgICAgbW9kZWwgPC0gdXBkYXRlKG1vZGVsLCBkYXRhPXRyKQogICAgICBwcmVkIDwtIHByZWRpY3QobW9kZWwsIG5ld2RhdGE9dmEpCiAgICAgIGlmIChuZWVkLmNsYXNzKXsKICAgICAgICBwcmVkIDwtIHByZWQkY2xhc3MKICAgICAgfQogICAgICAKICAgICAgZXJyLnRhYmxlIDwtIHRhYmxlKFRydWU9dGFyZ2V0W3ZhbF0sIFByZWRpY3RlZD1wcmVkKQogICAgICBlcnIgPC0gMS1zdW0oZGlhZyhlcnIudGFibGUpKS9zdW0oZXJyLnRhYmxlKQogICAgICBlcnIub25ldGltZSA8LSBjKGVyci5vbmV0aW1lLCBlcnIpCiAgICB9CiAgICBlcnIudG90YWwgPC0gYyhlcnIudG90YWwsIG1lYW4oZXJyLm9uZXRpbWUpKQogICAgcHJpbnQocGFzdGUwKCJJdGVyYXRpb24gIiwgaSwgIiwgbWVhbiBlcnJvcjogIiwgbWVhbihlcnIub25ldGltZSkpKQogIH0KICBtZWFuKGVyci50b3RhbCkKfQpgYGAKCgpDcm9zcy12YWxpZGF0aW9uIGZvciBrLU5lYXJlc3QgTmVpZ2hib3VyCgpgYGB7cn0KY3Jvc3MudmFsaWRhdGlvbi5rbm4gPC0gZnVuY3Rpb24oZGF0YSwgdGFyZ2V0LCB0aW1lcywgbmZvbGRzLCBLKXsKICBzZXQuc2VlZCgwMjAyKQogIENWLmZvbGRzIDwtIGdlbmVyYXRlQ1ZSdW5zKHRhcmdldCwgbnRpbWVzPXRpbWVzLCBuZm9sZD1uZm9sZHMsIHN0cmF0aWZpZWQ9VFJVRSkKCiAgZXJyLnRvdGFsIDwtIGMoKQogIGZvciAoaSBpbiAxOnRpbWVzKXsKICAgIGVyci5vbmV0aW1lIDwtIGMoKQogICAgZm9yIChqIGluIDE6bmZvbGRzKXsKICAgICAgdmFsIDwtIHVubGlzdChDVi5mb2xkc1tbaV1dW1tqXV0pCiAgICAgIAogICAgICB0ciA8LSBkYXRhWy12YWwsXQogICAgICB2YSA8LSBkYXRhW3ZhbCxdCgogICAgICBwcmVkIDwtIGtubih0ciwgdmEsIHRhcmdldFstdmFsXSwgayA9IEspCgogICAgICBlcnIudGFibGUgPC0gdGFibGUoVHJ1ZT10YXJnZXRbdmFsXSwgUHJlZGljdGVkPXByZWQpCiAgICAgIGVyciA8LSAxLXN1bShkaWFnKGVyci50YWJsZSkpL3N1bShlcnIudGFibGUpCiAgICAgIGVyci5vbmV0aW1lIDwtIGMoZXJyLm9uZXRpbWUsIGVycikKICAgIH0KICAgIGVyci50b3RhbCA8LSBjKGVyci50b3RhbCwgbWVhbihlcnIub25ldGltZSkpCiAgfQogIG1lYW4oZXJyLnRvdGFsKQp9CmBgYAoKCkNyb3NzLXZhbGlkYXRpb24gTmFpdmUtQmF5ZXMKCmBgYHtyfQpjcm9zcy52YWxpZGF0aW9uLm5haXZlIDwtIGZ1bmN0aW9uKGRhdGEsIHRhcmdldCwgbW9kZWwsIHRpbWVzLCBuZm9sZHMpewogIHNldC5zZWVkKDAyMDIpCiAgQ1YuZm9sZHMgPC0gZ2VuZXJhdGVDVlJ1bnModGFyZ2V0LCBudGltZXM9dGltZXMsIG5mb2xkPW5mb2xkcywgc3RyYXRpZmllZD1UUlVFKQoKICBlcnIudG90YWwgPC0gYygpCiAgZm9yIChpIGluIDE6dGltZXMpewogICAgZXJyLm9uZXRpbWUgPC0gYygpCiAgICBmb3IgKGogaW4gMTpuZm9sZHMpewogICAgICBwcmludChwYXN0ZTAoIkZvbGQ6ICIsIGopKQogICAgICB2YWwgPC0gdW5saXN0KENWLmZvbGRzW1tpXV1bW2pdXSkKICAgICAgCiAgICAgIHRyIDwtIGRhdGFbLXZhbCxdCiAgICAgIHZhIDwtIGRhdGFbdmFsLF0KCiAgICAgIG1vZGVsIDwtIG5haXZlQmF5ZXMoVjU4fi4tTEQxLUxEMi1MRDMtTEQ0LCBkYXRhPXRyKQogICAgICBwcmVkIDwtIHByZWRpY3QobW9kZWwsIG5ld2RhdGE9dmEpCiAgICAgIAogICAgICBlcnIudGFibGUgPC0gdGFibGUoVHJ1ZT10YXJnZXRbdmFsXSwgUHJlZGljdGVkPXByZWQpCiAgICAgIGVyciA8LSAxLXN1bShkaWFnKGVyci50YWJsZSkpL3N1bShlcnIudGFibGUpCiAgICAgIGVyci5vbmV0aW1lIDwtIGMoZXJyLm9uZXRpbWUsIGVycikKICAgIH0KICAgIGVyci50b3RhbCA8LSBjKGVyci50b3RhbCwgbWVhbihlcnIub25ldGltZSkpCiAgICBwcmludChwYXN0ZTAoIkl0ZXJhdGlvbiAiLCBpLCAiLCBtZWFuIGVycm9yOiAiLCBtZWFuKGVyci5vbmV0aW1lKSkpCiAgfQogIG1lYW4oZXJyLnRvdGFsKQp9CmBgYAoKCiMjIyA0LiBNb2RlbHMKVGhlIG1vZGVscyB3ZSBhcmUgZ29pbmcgdG8gdXNlIGFyZTogCiAgLSBMREEKICAtIFFEQQogIC0gUkRBCiAgLSBrLU5OCiAgLSBOYcOvdmUgQmF5ZXMKICAtIEdMTQogCgpgYGB7cn0KcmRhLm1vZGVsIDwtIHJkYShWNTh+VjMrVjQrVjkrVjEwK1YxMStWMTIrVjE0K1YxNStWMTYrVjE4K1YxOStWMjArVjIxK1YyMitWMjMrVjI0K1YyNStWMjYrVjI3K1YyOStWMzErVjMyK1YzMytWMzQrVjM1K1YzNytWMzgrVjM5K1Y0MCtWNDErVjQzK1Y0NCtWNTErVjU2K1Y2MCtWNjErVjYzK1Y2NStWNjcrVjY4LCBkYXRhPXRyYWluKQpuYWl2ZS5tb2RlbCA8LSBuYWl2ZUJheWVzKFY1OH5WMytWNCtWOStWMTArVjExK1YxMitWMTQrVjE1K1YxNitWMTgrVjE5K1YyMCtWMjErVjIyK1YyMytWMjQrVjI1K1YyNitWMjcrVjI5K1YzMStWMzIrVjMzK1YzNCtWMzUrVjM3K1YzOCtWMzkrVjQwK1Y0MStWNDMrVjQ0K1Y1MStWNTYrVjYwK1Y2MStWNjMrVjY1K1Y2NytWNjgsIGRhdGE9dHJhaW4pCmBgYAoKCmBgYHtyfQpjcm9zcy52YWxpZGF0aW9uKHRyYWluLCB0cmFpbiRWNTgsIHJkYS5tb2RlbCwgMTAsIDEwLCBUKQpgYGAKCgpgYGB7cn0KcmRhLm1vZGVsLmZkYSA8LSByZGEoVjU4fi4sZGF0YT10cmFpbikKYGBgCgoKYGBge3J9CmNyb3NzLnZhbGlkYXRpb24odHJhaW4sIHRyYWluJFY1OCwgcmRhLm1vZGVsLmZkYSwgMTAsIDEwLCBUKQpgYGAKCgpgYGB7cn0KY3Jvc3MudmFsaWRhdGlvbi5uYWl2ZSh0cmFpbiwgdHJhaW4kVjU4LCBuYWl2ZS5tb2RlbCwgMTAsIDEwKQpgYGAKCgpgYGB7cn0KZXJyIDwtIGMoKQpmb3IgKGsgaW4gMToyMCl7CiAgZXJyIDwtIGMoZXJyLCBjcm9zcy52YWxpZGF0aW9uLmtubih0cmFpbiwgdHJhaW4kVjU4LCAxMCwxMCwgaykpCn0KYGBgCgoKYGBge3J9CnBsb3QoZXJyLCB0eXBlID0gImwiKQplcnIKYGBgCgoKYGBge3J9CmNyb3NzLnZhbGlkYXRpb24ua25uKHRyYWluLCB0cmFpbiRWNTgsIDEwLCAxMCwgMSkKCm11bHRpbm9taWFsLm1vZGVsIDwtIG11bHRpbm9tKFY1OH4uLCBkYXRhPXRyYWluKQoKY3Jvc3MudmFsaWRhdGlvbih0cmFpbiwgdHJhaW4kVjU4LCBtdWx0aW5vbWlhbC5tb2RlbCwgMTAsIDEwLCBGKQoKbXVsdGlub21pYWwubW9kZWwuc3RlcCA8LSBzdGVwKG11bHRpbm9taWFsLm1vZGVsKQoKY3Jvc3MudmFsaWRhdGlvbih0cmFpbiwgdHJhaW4kVjU4LCBtdWx0aW5vbWlhbC5tb2RlbC5zdGVwLCAxMCwgMTAsIEYpCgptdWx0aW5vbWlhbC5tb2RlbC5ub0ZEQSA8LSBtdWx0aW5vbShWNTh+Li1MRDEtTEQyLUxEMy1MRDQsIGRhdGE9dHJhaW4pCgpjcm9zcy52YWxpZGF0aW9uKHRyYWluLCB0cmFpbiRWNTgsIG11bHRpbm9taWFsLm1vZGVsLm5vRkRBLCAxMCwgMTAsIEYpCgptdWx0aW5vbWlhbC5tb2RlbC5ub0ZEQS5zdGVwIDwtIHN0ZXAobXVsdGlub21pYWwubW9kZWwubm9GREEpCgpjcm9zcy52YWxpZGF0aW9uKHRyYWluLCB0cmFpbiRWNTgsIG11bHRpbm9taWFsLm1vZGVsLm5vRkRBLnN0ZXAsIDEwLCAxMCwgRikKYGBgCgoKIyMgVGVzdCBlcnJvcgoKYGBge3J9CnJkYS5tb2RlbCA8LSB1cGRhdGUocmRhLm1vZGVsLmZkYSwgZGF0YT10cmFpbikKcHJlZC50ZXN0IDwtIHByZWRpY3QocmRhLm1vZGVsLmZkYSwgdGVzdCkKcHJlZC50ZXN0IDwtIHByZWQudGVzdCRjbGFzcwooZXJyLnRhYmxlIDwtIHRhYmxlKFRydWU9dGVzdCRWNTgsIFByZWQ9cHJlZC50ZXN0KSkKKGVyci50ZXN0IDwtIDEtc3VtKGRpYWcoZXJyLnRhYmxlKSkvc3VtKGVyci50YWJsZSkpCmBgYAoKIyBQYXJ0ZSAyCgpXZSB3aWxsIG5vdyB1c2UgbW9kZWwgb24gYSBkYXRhc2V0IG5vdCB5ZXQgaW52ZXN0aWdhdGVkOiBIdW5nYXJ5LgoKYGBge3J9Cmh1bmcgPC0gcmVhZC50YWJsZSgiLi4vZGF0YS9wcm9jZXNzZWQuaHVuZ2FyaWFuLmRhdGEiLCBoZWFkZXI9Riwgc2VwPScsJywgbmEuc3RyaW5ncz0iPyIpCnN1bW1hcnkoaHVuZykKYGBgCgojIyBNaXNzaW5nIHZhbHVlcyB0cmVhdG1lbnQKCmBgYHtyfQpybS52YXIgPC0gYygiVjExIiwgIlYxMiIsICJWMTMiKQpodW5nIDwtIGh1bmdbLCEobmFtZXMoaHVuZykgJWluJSBybS52YXIpXQpgYGAKCmBgYHtyfQpodW5nW2lzLm5hKGh1bmcpXSA8LSAtOQp2YXIuaW1wdXRhYmxlIDwtIGMoIlY0IiwgIlY1IiwgIlY2IiwgIlY3IiwgIlY4IiwgIlY5IikKZm9yIChub20gaW4gdmFyLmltcHV0YWJsZSl7CiAgaHVuZ1ssIG5vbV1baHVuZ1ssIG5vbV0gPT0gLTldIDwtIE5BCiAgdmFyaWFibGUgPC0gaHVuZ1ssIG5vbV0KICBodW5nWywgbm9tXSA8LSBrbm4uaW1wdXRhdGlvbihodW5nLCB2YXJpYWJsZSwgbm9tLCA3KQp9CmBgYAoKYGBge3J9CnZhci5mYWN0b3IgPC0gYygiVjIiLCAiVjMiLCAiVjYiLCAiVjciLCAiVjkiLCAiVjE0IikKZm9yIChub20gaW4gdmFyLmZhY3Rvcil7CiAgaHVuZ1ssIG5vbV0gPC0gYXMuZmFjdG9yKGFzLmNoYXJhY3RlcihodW5nWyxub21dKSkKfQoKIyBMZXZlbCAyIGhhcyBmZXcgdmFsdWVzCmh1bmckVjdbaHVuZyRWNyA9PSAyXSA8LSAxCmh1bmckVjcgPC0gZHJvcGxldmVscyhodW5nJFY3KQpsZXZlbHMoaHVuZyRWNykKCiMgTGV2ZWwgMiBoYXMgb25lIGFub21hbCB2YWx1ZQpodW5nJFY5W2h1bmckVjkgPT0gMl0gPC0gMQpodW5nJFY5IDwtIGRyb3BsZXZlbHMoaHVuZyRWOSkKbGV2ZWxzKGh1bmckVjkpCmBgYAoKCmBgYHtyfQpzdW1tYXJ5KGh1bmcpCmBgYAoKIyMgR2F1c3NpYW5pdHkgdHJlYXRtZW50CgpgYGB7cn0KbmFtZXNfbnVtIDwtIGMoKQpmb3IoaSBpbiAxOm5jb2woaHVuZykpewogIGlmICghaXMuZmFjdG9yKGh1bmdbLGldKSkgewogICAgbmFtZXNfbnVtIDwtIGMobmFtZXNfbnVtLCBpKQogIH0KfQpodW5nX251bWVyaWMgPC0gaHVuZ1ssbmFtZXNfbnVtXQoKaHVuZ19jb3IgPC0gY29yKGh1bmdfbnVtZXJpYykKd2hpY2goaHVuZ19jb3IgPiAwLjUgJiBodW5nX2NvciA8IDEsIGFyci5pbmQgPSBUUlVFKQp3aGljaCgtaHVuZ19jb3IgPiAwLjUgJiAtaHVuZ19jb3IgPCAxLCBhcnIuaW5kID0gVFJVRSkKYGBgCgoKYGBge3J9CnBhcihtZnJvdyA9IGMoMiwzKSkKZm9yKGkgaW4gMTpsZW5ndGgoaHVuZ19udW1lcmljKSl7CiAgICBxcW5vcm0oaHVuZ19udW1lcmljWyxpXSwgbWFpbiA9IGMoIlEtUSBQbG90OiAiLCBuYW1lcyhodW5nX251bWVyaWMpW2ldKSkKICAgIHFxbGluZShodW5nX251bWVyaWNbLGldLCBjb2w9MikKfQpgYGAKCgpCb3hjb3gKCmBgYHtyfQpwYXIobWZyb3cgPSBjKDIsMykpCmZvcihpIGluIDE6bGVuZ3RoKGh1bmdfbnVtZXJpYykpewogIGJveGNveChsbShodW5nX251bWVyaWNbLGldLW1pbihodW5nX251bWVyaWNbLGldKSsxZS02fjEpLGxhbWJkYSA9IHNlcSgtMSwgMS41LCBieT0wLjEpLCB4bGFiID0gYyhuYW1lcyhodW5nX251bWVyaWMpKVtpXSkKfQpgYGAKCldlIHNlZSB2YXIgMTAgcmVxdWlyZXMgYSBsb2dhcml0bWljIHRyYW5zZm9ybWF0aW9uLgoKYGBge3J9Cmh1bmdbLCJWMTAiXSA8LSBsb2coMStodW5nWywiVjEwIl0pCmBgYCAKCmBgYHtyfQpxcW5vcm0oaHVuZ1ssIlYxMCJdKQpxcWxpbmUoaHVuZ1ssIlYxMCJdLCBjb2w9MikKYGBgCgojIyBGZWF0dXJlIGV4dHJhY3Rpb24KClNlcGFyZSB0cmFpbiBhbmQgdGVzdCBkYXRhLCBzZWVkIGZvciByZXByb2R1Y2liaWxpdHkuCgpgYGB7cn0Kc2V0LnNlZWQoMjAwMCkKbiA8LSBucm93KGh1bmcpCnRyYWluLmxlbmdodCA8LSByb3VuZCgyKm4vMykKCmh1bmcgPC0gaHVuZ1tzYW1wbGUobiksXQp0cmFpbiA8LSBodW5nWzE6dHJhaW4ubGVuZ2h0LF0KdGVzdCA8LSBodW5nWyh0cmFpbi5sZW5naHQrMSk6bixdCmBgYAoKYGBge3J9Cm5hbWVzX251bSA8LSBjKCkKZm9yKGkgaW4gMTpuY29sKHRyYWluKSl7CiAgaWYgKCFpcy5mYWN0b3IodHJhaW5bLGldKSkgewogICAgbmFtZXNfbnVtIDwtIGMobmFtZXNfbnVtLCBpKQogIH0KfQp0cmFpbl9udW0gPC0gdHJhaW5bLG5hbWVzX251bV0KYGBgCgpFeHRyYWN0IFBDQSBmZWF0dXJlcy4KCmBgYHtyfQpwY2EgPC0gcHJpbmNvbXAodHJhaW5fbnVtKQpzY3JlZXBsb3QocGNhKQpzdW1tYXJ5KHBjYSkKYGBgCgoKYGBge3J9CmJpcGxvdChwY2EpCmBgYAoKYGBge3J9CkZwIDwtIHBjYSRzY29yZXMKR3MgPC0gcGNhJGxvYWRpbmdzCgpGcyA8LSBGcCAlKiUgZGlhZygxL3BjYSRzZGV2KQpHcCA8LSBHcyAlKiUgZGlhZyhwY2Ekc2RldikgKjAuMwoKY29sLmNsYXNzIDwtIGFzLm51bWVyaWModHJhaW4kVjE0KQpjb2wuY2xhc3NbY29sLmNsYXNzPT0wXSA8LSAicmVkIgpjb2wuY2xhc3NbY29sLmNsYXNzPT0xXSA8LSAiYmx1ZSIKCnBsb3QoRnNbLDFdLCBGc1ssMl0sIGFzcD0xLCBjb2wgPSBjb2wuY2xhc3MsIHhsYWIgPSAiRmlyc3QgcHJpbmNpcGFsIGNvbXBvbmVudCIsIHlsYWIgPSAiU2Vjb25kIHByaW5jaXBhbCBjb21wb25lbnQiKQphcnJvd3MocmVwKDAsZGltKEdzKVsxXSkscmVwKDAsZGltKEdzKVsxXSksIEdwWywxXSwgR3BbLDJdKQp0ZXh0KEdwWywxXSwgR3BbLDJdLCBuYW1lcyh0cmFpbl9udW0pLCBjb2wgPSAiYmxhY2siKQpsZWdlbmQoImJvdHRvbXJpZ2h0IiwgZmlsbD1jKCJyZWQiLCJibHVlIiksIGxlZ2VuZD1jKCcwJywnMScpKQpgYGAKCmBgYHtyfQpmZGEgPC0gbGRhKFYxNH4uLCBkYXRhPXRyYWluKQojcGxvdChmZGEpCmxvYWRpbmdzIDwtIHByZWRpY3QoZmRhKSR4CnBsb3QobG9hZGluZ3MsIGNvbCA9IGNvbC5jbGFzcykKbGVnZW5kKCJib3R0b21yaWdodCIsIGZpbGw9YygicmVkIiwiYmx1ZSIpLCBsZWdlbmQ9YygnMCcsJzEnKSkKYGBgCgoKYGBge3J9CnRyYWluJExEMSA8LSBsb2FkaW5nc1ssMV0KCmZkYV90ZXN0IDwtIHByZWRpY3QoZmRhLCBuZXdkYXRhID0gdGVzdCkKdGVzdCRMRDEgPC0gZmRhX3Rlc3QkeFssMV0KYGBgCgoKQW7DoWxpc2lzIHBvciBjb3JyZXNwb25kZW5jaWFzCgpgYGB7cn0KYWMgPC0gbWpjYShodW5nWyxuYW1lcyhodW5nKSAlaW4lIHZhci5mYWN0b3JdLCBsYW1iZGE9IkJ1cnQiKQpwbG90KGFjLCBtYWluPSJNQ0EgYmlwbG90IG9mIEJ1cnQgbWF0cml4IHdpdGggZGF0YSIpCgphY19pbmQgPC0gbWpjYSh0cmFpblssbmFtZXModHJhaW4pICVpbiUgdmFyLmZhY3Rvcl0sIGxhbWJkYT0iaW5kaWNhdG9yIiwgcmV0aSA9IFQpCnBsb3QoYWNfaW5kJHJvd2Nvb3JkLCBjb2wgPSBjb2wuY2xhc3MpCmxlZ2VuZCgiYm90dG9tcmlnaHQiLCBmaWxsPWMoInJlZCIsImJsdWUiKSwgbGVnZW5kPWMoJzAnLCcxJykpCgptY2EuZmVhdHVyZXMgPC0gYWNfaW5kJHJvd2Nvb3JkCmBgYAoKIyMgCg==